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Abstract 

An  efficient  numerical  scheme  is  presented  for  simulating  the 
isothermal  flow  in  resin  transfer  molding  (RTM).  The  problem  involves 
transient,  free  surface  flow  of  an  incompressible  fluid  into  a  nonde¬ 
forming  porous  medium.  A  new  variant  of  the  control  volume  finite 
element  (CVFE)  algorithm  is  explained  in  detail.  It  is  shown  how  the 
pressure  solutions  at  each  time  step  can  be  obtained  by  adding  a  single 
row  and  column  to  the  Cholesky  factorization  of  the  stiffness  matrix 
derived  from  a  finite-element  formulation  for  the  pressure  field.  This 
approach  reduces  the  computation  of  a  new  pressure  solution  at  each 
time  step  to  essentially  just  two  sparse  matrix  back-substitutions.  The 
resulting  performance  improvement  facilitates  interactive  simulation 
and  the  solution  of  inverse  problems  which  require  many  simulations 
of  the  filling  problem.  The  computational  complexity  of  the  calcula¬ 
tion  is  bounded  by  where  n  is  the  number  of  nodes  in  the 

finite  element  mesh.  A  100-fold  speedup  over  a  conventional  CVFE 
implementation  was  obtained  for  a  2213-node  problem. 
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1  Resin- Transfer  Mold  Filling  Problem 

Resin  transfer  molding  (RTM)  is  an  emerging  manufacturing  technology  well- 
suited  for  fabricating  large  structural  components  made  of  composite  mate¬ 
rials.  Since  the  process  involves  matched  metal  tooling,  the  technique  seems 
ideal  for  situations  requiring  close  tolerances.  Construction  of  aircraft  struc¬ 
tures  and  vehicle  components  fit  this  characterization.  Furthermore,  liquid 
injection  molding  represents  one  of  the  most  economical  means  of  manufac¬ 
turing.  RTM  is  an  adaptation  on  a  process  widely  used  for  plastics.  Instead 
of  injecting  into  an  empty  cavity,  the  mold  is  packed  with  a  woven  fiber 
preform.  The  RTM  process  has  two  main  stages:  filling  the  mold  with  a 
resin/catalyst  mixture  and  curing  the  part. 

At  present,  most  of  the  difficulties  of  incorporating  RTM  revolving  around 
the  filling.  To  create  an  acceptable  composite  part  requires  the  preform  to 
be  completely  impregnated  with  resin.  This  largely  controlled  by  the  fluid 
dynamics  of  the  resin  flow  into  the  fiber  reinforcement.  The  conditions  which 
most  strongly  influence  the  the  flow  are:  mold  geometry,  resin  rheology,  pre¬ 
form  permeability,  and  location  of  the  injectors  and  vents.  The  first  three 
are  typically  determined  by  the  part  design  itself;  the  last  one  is  a  manu¬ 
facturing  consideration.  Incorrect  placement  of  the  injectors  and  vents  for  a 
given  geometry  and  resin/preform  system  will  create  dry  spots  in  the  cured 
part. 

To  enhance  the  economic  viability  of  RTM  applications,  it  is  desirable  to 
evaluate  mold  design  prior  to  mold  construction,  via  computer-based  meth¬ 
ods.  Predictive  modeling  of  resin  flow  through  a  fiber  preform  is  currently 
an  important  priority  in  mold  design  and  evaluation,  because  of  the  need 
to  predict  fill  times  and  wet-out  patterns.  Ongoing  research  also  includes 
control  algorithms  for  the  filling  stage,  using  networks  of  embedded  sensors 
and  a  fast  filling  simulation  [1,2].  ’ 

Simulation  of  the  RTM  process  may  be  isothermal  or  nonisothermal,  de¬ 
pending  upon  whether  temperature  effects  are  accounted  for  in  the  model. 
During  filling,  resin  viscosity  is  affected  by  temperature  variations.  During 
curing,  gel  times  are  affected  by  temperature  profiles.  For  a  mathematical 
formulation  of  the  nonisothermal  RTM  process,  see  [4].  We  focus  on  the 
isothermal  case  under  the  assumption  of  minimal  temperature  variation  dur¬ 
ing  filling. 

The  isothermal  RTM  filling  problem  is  a  transient,  free-boundary  prob- 
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lem  of  predicting  the  position  of  the  resin  flow  front  in  the  porous  medium 
as  a  function  of  injection  pressures  and  time.  The  resin  is  assumed  to  be 
nearly  compressible  and  to  display  Newtonian  behavior.  The  fiber  preform 
is  assumed  to  be  non- deforming.  It  is  assumed  that  Darcy’s  law  governs  the 
relation  between  resin  velocity,  v,  and  pressure  p,  such  that 

v  =  -p-^KVp,  (1) 

where  fi  is  the  viscosity  and  K  is  a  tensor  representing  the  permeability  of 
the  fiber  preform.  Preforms  are  usually  constructed  from  several  layers  of 
fiber  mat  oriented  in  different  directions.  Permeabilities  are  experimentally 
calculated  for  mat  samples  and  reported  in  terms  of  the  principal  directions 
of  the  mat.  Thus,  the  average  (through-thickness)  permeability  is  a  function 
of  several  factors,  see  for  example  [8].  Since  permeabilities  along  the  principal 
axes  can  easily  differ  by  an  order  of  magnitude,  the  ability  to  specify  K  on 
a  local  basis  is  essential  in  simulation. 

Let:  12  C  define  the  interior  of  the  mold;  T^u  the  impermeable  mold 
walls;  P/i  the  constant  displacement  injectors;  and  P^  the  constant  pressure 
injectors;  so  that  the  complete  assembly  includes 

12  =  i2ur„ur,,urg.  (2) 

Let:  fl{t)  C  12  denote  the  filled  portion  of  the  mold  interior  at  time  t  and 
Ts{t)  the  free  surface  at  time  t.  On  the  air  side  of  the  surface,  the  capillary 
fringe  is  neglected  and  a  constant  pressure  (e.g.  atmospheric)  is  assumed  in 
the  unfilled  portion  of  the  mold  {12  \  f2(f)}. 

The  isothermal  RTM  filling  problem  is  to  find  for  any  t  >  0,  Ts{t)  :  h-> 

and  p  :  12(2)  i— R  such  that 

V-(^-i/PVp) 
n  •  {p-^KVp) 
n  •  {p-^KVp) 

P 
P 

n  ■  (p-^KVp) 


0 

in  12(2) 

(3) 

0 

on 

(4) 

h{x) 

on  Vh 

(5) 

gi^) 

on  Tg 

(6) 

0 

on  rs(2) 

(7) 

dTs 
— n  •  — — 
dt 

.  -  . 

(8) 
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where  n  denotes  the  vector  normal  to  F.  The  quantity  /  =  n  •  v  is  defined  as 
the  flux  normal  to  the  boundary.  The  mold-  filling  problem  is  analogous  to 
the  Stefan  problem,  a  class  of  free  boundary  problems  used  in  modeling  the 
melting  of  solids  and  crystallization  of  liquids  [3,  15]. 

2  Filling  Algorithm 

The  mold  filling  problem  in  equations  (3)  -  (8)  may  be  solved  with  a  variety 
of  numerical  schemes,  including  fully  implicit  methods  which  solve  simulta¬ 
neously  for  free  surface  location  and  pressure,  as  well  as  more  conventional 
semi-implicit  methods  which  solve  for  pressure  implicitly  and  satisfy  the  free 
surface  condition  (8)  by  an  explicit  method.  We  follow  the  latter  approach, 
suggested  by  many  researchers  [5,  9,  4,  6],  using  a  finite  element  solution  to 
obtain  pressures  at  discrete  values  of  t  followed  by  an  explicit  control  volume 
scheme  for  updating  the  position  of  the  free  surface.  The  control  volume 
scheme  involves  domain  discretization  into  discrete  subvolumes,  where  each 
such  control  volume  contains  one  node  of  the  finite  element  mesh.  Fluid  flux 
across  control  volume  boundaries  is  calculated  from  the  pressure  solution. 
Net  inflow  to  a  control  volume  is  tracked  as  the  “fill  fraction”.  When  the 
volume  is  filled,  the  node  contained  in  the  volume  is  considered  part  of  the 
next  pressure  solution,  as  summarized  below. 

1.  Find  FEM  pressure  solution.  At  the  beginning  of  a  new  time  step,  the 
pressure  field  is  calculated  over  the  filled  control  volumes. 

2.  Calculate  volume  flux.  Darcy  velocity  and  flux  at  the  control  volume 
boundaries  are  computed  from  the  pressures  calculated  in  step  1. 

3.  Locate  the  free  surface.  The  time  step  is  calculated  as  the  minimum 
At  required  to  fill  a  control  volume  using  the  flux  calculated  in  step  2. 
The  filled  volume  becomes  part  of  the  pTessure  field,  moving  the  free 
surface  to  a  new  position. 

2.1  Finite  Element  Pressure  Solution 

At  each  time  step,  a  free  surface  location  Ts{t)  is  given  from  the  .previous  step. 
The  problem  of  solving  equation  (3)  subject  to  (4)-(7)  is  analogous  to  the 
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classical  heat  conduction  problem,  which  has  the  following  weak  formulation. 
Define  the  trial  function  and  weight  function  spaces. 


U  =  {p  I  p  e p|r,(f)  =  p|r,  =  ^(x)} 
V  =  {w  \  w  e  II^,  u^lrjur,  =  0). 

Given  g,  h,  find  p  G  U  such  that  for  all  w  G  V, 


[  (Vwf(g-'^KWp)dn  =  [  whdTh. 
Jq 


The  Galerkin  formulation  (in  two  spatial  dimensions)  is  based  on  a  tri¬ 
angulation  of  the  mold  into  nel  elements, 

f2  «  T  =  UiTi,  i  =  1,  ...,nel, 

where  the  Ti  are  closed  triangles  with  nodes  (vertices)  Xj  =  {xj,yj),j  — 
(see  Figure  2). 

We  have  chosen  to  work  in  the  space  of  piecewise  linear  functions.  Let  p 
be  the  finite  dimensional  approximation  to  p,  such  that  p  G  U  is  &  polynomial 
of  degree  one  over  each  triangle.  The  vector  p  will  denote  the  nodal  values 
of  p,  i.e., 


[Pi,-,P. 


(Xi),...,p(Xn) 


and  p(x)  =  A^i(x)pi),  where  the  shape  functions  Ni{xj)  = 
1,  ...,n  are  piecewise  linear.  The  Galerkin  formulation  is  then 


J2(l  ViV,  •  (p-^/FViV,)dD)  p,  =  [  NMTh,  z  =  l,...,n. 

\Jq  J  JVh  ' 

The  elements  of  the  n  x  n  pressure  stiffness  matrix  A  are  defined 

aij  =  f  VNi  ■  {p~^K'VNj)dfl,  i,j  = 

Jq 

and  the  forcing  vector  b  as 

bi  =  /  NihdVf,,  i  =  l,...,n. 

Jr^ 
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In  matrix  notation. 


Ap  =  h,  (14) 

and  the  approximate  pressure  solution  is  obtained  by  solving  a  linear  system 
of  n  equations  in  n  unknowns. 

2.2  Control  Volume  Flux  Calculation 

The  flux  calculation  approximates  fluid  velocities  by  the  piecewise  linear 
polynomials  described  in  the  previous  section, 

v(x)  =  X]/^''AV7V,(x)pi. 

!=1 

We  use  a  node-centered  control  volume  discretization  to  calculate  flux 
and  fluid  volume  fractions.  The  control  volume  discretization  is  built  on  the 
triangulation  T  (see  Figure  2)  and  forms  a  set  of  closed  subvolumes  Ci, 

ft  Ri  C  =  UjCj,  i  =  1, ...,  n, 

such  that  subvolume  Ci  contains  vertex  Xj-  and  no  other  node.  Denoting  Bi 
as  the  boundary  of  C,-,  the  flow  rate  into  Ci  at  time  tis 

qi{t)  =  /  n  ■  v{t)dBi.  (15) 

J  Bi 

where  n  is  the  unit  normal  to  Bi. 

We  define  Bi  as  the  set  of  line  segments  connecting  element  centroids  with 
edge  midpoints,  as  follows  (see  Figure  3).  Associated  with  Bi  is  the  set  of 
elements  Ei  =  {Tj  \  x,  G  Tj}.  Without  loss  of  generality,  denote  the  vertices 
of  Tj  as  X,-, Xj+i, Xi+2.  Define  the  centroid  of  Tj  as  cj  =  \  and  the 

edge  midpoints,  myi  =  |(x, +x,+i)  and  mj2  =  |(xi  +  x,+2)-  The  segments 
of  Bi  in  Tj  are  denoted  by  (cj,mji),  (cj,mj2).  Then  qi  can  be  written  as  an 
algebraic  sum. 
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*li  —  ^  (/jinji  -I-  Ij2nj2)  ■  ^  fi^^KVNkPk  (16) 

TjEEx  k=l 

where  n^i  is  the  unit  normal  vector  in  the  plane  of  Tj  orthogonal  to  (cj,mji), 
and  Iji  =  ||cj  —  mji||2  is  the  segment  length.  The  coefficients  of-p  in  equation 
(16)  are  constant  for  all  t  and  are  assembled  prior  to  the  filling  algorithm. 
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Figure  2:  Triangulation  and  Control  Volume  Discretization  of  O 
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■  Figure  3:  Control  Volume  Discretization,  (a)  The  set  Ei  of  elements 
supported  by  node  i,  (b)  centroid  and  edge  .midpoints  for  Ti,  (c)  control 
volume  Ci  and  bouundary  B,,  (d)  superimposed  FEM  and  control  volume 
discretizations. 
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2.3  Free  Surface  Location 

A  node  X,-  is  included  in  the  filled  domain  Cl{t)  if  control  volume  Ci  has  a  fill 
fraction  of  unity.  Let  Si{t)  denote  the  fill  fraction  of  Ci  at  time  f, 

Si{t)  =  \Ci\-^  (17) 

Jo 

where  \Ci\  denotes  volume  adjusted  for  the  porosity  of  the  preform.  At  each 
time  step,  the  fill  fraction  is  updated  explicitly.  If  qi{t)  is  the  flow  rate  into 
Ci  at  time  f,  then 


Si{t  +  M)  =  S,{t)  +  Mq,{t).  (18) 

According  to.  equation  (7),  the  free  surface  location  satisfies  the  relation 
pItj  =  0.  We  define  the  free  surface  location  by  the  nodes  in  unfilled  or 
partially  filled  volumes  adjacent  to  filled  nodes,  i.e.,  Xj  is  on  the  free  surface 
if  (a)  Sj  <  1,  and  (b)  Xj  is  adjacent  to  a  node  x^  such  that  =  1.  Observe 
in  Figure  2  the  free  surface  intersects  partially  filled  control  volumes  and 
that  control  volume  flux  is  not  actually  calculated  at  the  free  surface.  To 
address  this  fundamental  discrepancy,  some  researchers  have  developed  local 
refinement  schemes  in  the  flow  front  vicinity,  e.g.  [10].  Such  schemes  improve 
the  local  accuracy  of  the  flow  front  approximation.  However,  we  believe  that 
if  high  accuracy  in  flow  front  calculations  is  needed,  then  an  alternative  filling 
algorithm  satisfying  more  rigorous  mathematical  convergence  criteria  should 
be  considered  rather  than  local  mesh  refinement. 

2.4  Implementation 

2.4.1  Properties  of  the  Pressure  Stiffness  Matrix 

The  pressure  stiffness  matrix  A  in  equation  (14)  has  several  important  prop¬ 
erties  which  lead  to  an  efficient  implementation  of  the  filling  algorithm.  The 
first  property  is  sparsity,  due  to  the  structure  of  the  finite  element  mesh. 
The  order  of  A  is  equal  to  the  number  of  nodes,  n  in  the  mesh.  The  number 
of  nonzero  entries  in  A  is  equal  to  the  number  of  edges  connecting  adjacent 
nodes.  Since  the  nodes  of  a  finite  element  mesh  typically  are  connected  to 
only  a  few  other  nodes,  the  number  of  nonzeroes  is  far  less  than  n^,  usually 
a  small  multiple  of  n.  The  second  and  third  properties  of  A  are  symmetry, 
A  =  A^  and  positive-definiteness,  A(A)  >  0. 
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To  demonstrate  positive-definiteness,  it  is  necessary  to  show  that  K  is 
always  positive-definite.  Experimental  permeability  measurements  are  re¬ 
ported  for  principal  mat  directions  as  a  diagonal  matrix  with  strictly  positive 
entries,  D  =  diag{k\\^k22ikzz).  The  tensor  K  is  obtained  by  transforming 
from  the  principal  directions  of  the  mat  to  the  Cartesian  frame  (or  to  local 
element  coordinates),  i.e.,  K  =  CDC^,  where  (7  is  a  rank-3,  orthogonal  ro¬ 
tation  matrix  which  projects  the  principal  axes  of  the  mat  into  the  Cartesian 
or  local  coordinate  system.  Since  pre  or  post  multiplication  by  an  orthogo¬ 
nal  matrix  preserves  the  spectrum  of  an  operator,  A(A")  =  X{D)  >  0.  Given 
A(A^)  >  0,  symmetry  and  positive  definiteness  of  A  is  a  standard  result  in 
the  finite  element  literature;  see,  for  example,  [14]. 

Since  A  is  symmetric  and  positive  definite,  it  has  a  Cholesky  factorization, 
LL^  =  A,  where  A  is  a  lower  triangular  matrix.  It  is  also  true  that  every 
submatrix  of  A  inherits  these  two  properties.  Thus,  if  A  is  partitioned  as 

.  (  M  u\ 

^  s  j 

then  the  matrix  M  has  a  Cholesky  factorization  M  =  LmL'm 

where  Lm'^^  =  u  and  t  =  [s  —  vo^w)^ .  As  a  result,  the  Cholesky  factor  of  the 
stiffness  matrix  can  be  computed  row  by  row.  This  is  exactly  the  property 
required  for  an  efficient  isothermal  filling  algorithm. 

2.4.2  Updating  the  Pressure  Solution 

Each  time  step  At  is  calculated  to  fill  one  control  volume.  Filled  volumes 
are  considered  part  of  the  fluid  phase  and  so  the  corresponding  node  be¬ 
comes  part  of  the  fluid  pressure  calculation.  The  addition  of  a  node  to  the 
pressure  calculation  corresponds  to  adding  a  single  row  and  column  to  the 
stiffness  matrix.  The  Cholesky  factorization  of  the  updated  stiffness  matrix 
can  be  updated  directly,  as  in  equation  (20).  The  advantage  is  that  the  stiff¬ 
ness  matrix  A  need  only  be  factored  one  time,  rather  than  reassembling  and 
factorizing  at  every  time  step. 
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The  stiffness  matrix  A  is  assembled  and  stored  prior  to  filling.  The  full 
matrix  is  stored  in  an  adjacency  structure.  The  adjacency  structure  consists 
of  n  adjacency  lists  and  corresponding  nonzero  coefficients.  The  ith  adja¬ 
cency  list  includes  the  indices  of  nodes  which  are  adjacent  to  (share  an  edge 
with)  node  x*-. 

During  the  filling  algorithm,  nodes  are  added  to  the  pressure  field  as 
control  volumes  are  filled.  Rows  of  A  corresponding  to  these  nodes  are  added 
to  the  Cholesky  factor  L  using  Equation  (20).  Note  that  these  rows  must  be 
permuted  to  reflect  the  node  ordering  imposed  by  the  filling  sequence.  The 
data  structure  used  to  store  L  is  an  envelope  structure.  For  each  row  of  the 
matrix,  all  entries  from  the  first  nonzero  up  to  the  diagonal  are  stored. 

The  envelope  storage  scheme  is  a  standard  data  structure  for  sparse  ma¬ 
trix  factorization.  This  choice  permits  the  use  of  existing  numerical  soft¬ 
ware  [12]  for  updating  the  Cholesky  factorization  and  computing  intermedi¬ 
ate  pressure  solutions  at  each  time  step,  with  only  minor  modifications. 

In  practice,  the  algorithm  will  often  fill  more  than  one  node  in  a  single 
time  step,  despite  the  fact  that  the  time  step  is  calculated  to  fill  only  one 
volume.  This  occurs  most  typically  in  regular  discretizations  because  a  tol¬ 
erance  is  used  to  define  the  fill  fraction  constituting  a  “filled”  control  volume 
(e.g.,  99%). 

2.4.3  CVFE  Algorithm 

The  CVFE  algorithm  requires  an  extensive  set  of  inputs  and  initialization 
steps.  For  the  case  of  a  two-dimensional  thin  shell  geometry  in  three  dimen¬ 
sional  space,  these  initialization  steps  include: 

•  specification  of  a  triangulation  T, 

•  specification  of  a  control  volume  discretization  C, 

•  specification  of  local  permeabilities  and” element  thickness, 

•  rotation  of  permeabilities  to  local  element  coordinates, 

•  calculation  of  adjacency  data  structure  for  A, 

•  assembly  of  A  in  adjacency  structure,  -  . 

•  assembly  of  flow  rate  coefficient  matrix. 
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The  iterative  part  of  the  algorithm  computes  the  filling  sequence,  using 
the  integer  arrays  perm  and  invp  to  denote  the  ordering  of  A  (the  filling 
sequence)  used  in  the  Cholesky  factorization.  The  notation  perm(i)  =  k 
means  the  original  node  is  the  zth  node  in  the  new  ordering.  The  element 
invp{k)  gives  the  position  in  perm  where  the  node  originally  numbered  k 
resides,  i.e.  perm{invp{k))  =  k.  We  use  the  array  subscript  notation  hi  =  h{i) 
interchangeably. 

^  rr  0;  /  =  0;  m  =  0;  filled  =  false; 
while  (not  filled) 

compute  forcing  vector  b 

{add  row(s)  to  pressure  stiffness  factor  L} 

for  z  =  1, ...,  n  {add  filled  control  volumes  to 

if  Si{t)  =  1  and  x,-  3  Q{t), 

m  =  m  +  1;  {increment  number  of  filled  volumes} 

perm{i)  =  m; 
invp{m)  =  z; 

{scatter  row  z  from  adjacency  structure  of  A  to  full  vector} 

{gather  permuted  row  i  into  envelope  structure  for  L} 

Lij  w{invp{j)),  j  =  l,...,n; 
end  if 
end  for 

{solve  updated  pressure  system} 

b(perm(z))  ^  b(z),  z  =  1, ...,  n;  {permute  forcing  vector} 
update  Xij,  z  =  /+l,...,m,  {e'qn  (20)} 

solve  y  iy  =  b; 
solve  p  ^  =  y; 

I  =  m;  {updated  dimension  of  L} 

p{invp{i))  ^  p^,  z  ==  1, ...,  m;  {restore  solution  to  original  order} 

{update  fill  fractions} 

compute  qi{t)^  i  =  1, ...,  n; 
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if  qi{t)  ==  0,  i  ='l,  STOP;  {mold  cannot  be  filled} 
At  =  mini[(l  -  Si{t))  \Ci\  /  i  =  l,...,n; 

Si{t  +  At)  =  Si{t)  +  Atqi{t),  i  =  1, 71] 
t  =  t  +  At] 

if  Si{t)  =  1,  i  =  1, n,  filled  ~  true; 
end  while 


The  order  in  which  control  volumes  are  filled  determines  the  ordering  of 
A  during  the  Cholesky  factorization.  A  more  efficient  ordering  could  be  ob¬ 
tained  by  using  a  symbolic  factorization  algorithm,  such  as  reverse  Cuthill- 
McKee  [11]  to  find  an  ordering  which  reduces  the  maximum  envelope  band¬ 
width.  Such  a  procedure  reduces  the  number  of  nonzero  elements  in  L  and 
hence  the  computational  effort  in  solving  Ly  =  b  and  L^x  =  y.  This  possi¬ 
bility  will  be  addressed  in  the  context  of  further  research  on  the  application 
of  direct  factorization  methods  for  large-scale  RTM  simulation  and  parallel 
computation. 

2.5  Computational  Complexity 

This  section  analyzes  the  computational  requirements  of  the  filling  algorithm 
as  a  function  of  the  problem  size.  The  analysis  is  based  on  several  assump¬ 
tions  about  typical  models  and  is  not  a  worst-case  analysis  of  complexity. 

The  CVFE  algorithm  requires  0{n)  iterations  or  time  steps,  one  per 
control  volume  (in  practice,  more  than  one  volume  may  fill  per  time  step). 
Each  iteration  requires  the  four  procedures  as  summarized  below. 

for  A;  =  1, ...,  n 

compute  forcing  vector  b 
add  row(s)  to  L 
solve  updated  pressure  system 
update  fiU  fractions 
end  while 


0{h)  operations 
-  0{n)  operations 

operations 
0{n)  operations 


Three  procedures  require  0{n)  operations  per  iteration,  or  0{n‘^)  opera¬ 
tions  for  all  iterations.  One  procedure,  solving  the  updated  pressure  system. 
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requires  Oik^'^)  operations  per  iteration,  where  k  is  the  iteration  number. 
This  procedure  dominates  the  total  computation  and  is  explained  below. 

The  updated  pressure  system  requires  several  sparse-matrix  back  substi¬ 
tutions  using  L.  The  number  of  operations  in  a  back  substitution  is  propor¬ 
tional  to  the  number  of  nonzeroes  in  L.  It  is  assumed  that  the  A  originally 
has  0{n)  nonzeroes.  This  is  reasonable  since  the  maximum  adjacency  list 
length  is  typically  a  small  constant  (e.g.,  6).  However,  it  is  also  assumed 
that  A  has  the  structure  of  an  n  x  n  Laplacian  matrix,  with  a  bandwidth 
of  0{\/n),  and  that  no  reordering  scheme  is  used.  Then  the  fill-in  of  L  is 
nonzeroes,  corresponding  to  fill-in  between  the  bands.  It  is  assumed 
that  fill-in  occurs  at  the  rate  0{k^'^)  per  iteration  as  rows  are  added  to  L. 
Using  the  following  relationship, 

k^-^  <y/^^k=  ^\/n(n^  -1-  n), 

k=l  k=l  ^ 

we  bound  the  total  computational  effort  as  0(n^-®).  Thus,  the  complete 
CVFE  algorithm  requires  operations.  A  CVFE  algorithm  which  as¬ 

sembles  and  factors  the  stiffness  matrix  at  each  iteration  would  theoretically 
require  O(n^)  operations  per  time  step  and  0{n^)  in  total.  The  introduc¬ 
tion  of  reordering  algorithms  could  further  reduce  the  complexity  of  both  the 
CVFE  algorithm  described  in  this  paper  as  well  as  conventional  approaches. 

3  Numerical  Results 

This  section  summarizes  implementation  and  performance  details  of  the  fill¬ 
ing  algorithm.  The  details  include  type  of  architecture  and  source  language, 
method  of  validation,  comparison  with  related  codes,  and  timing  results  for 
several  test  problems. 

The  filling  algorithm  has  been  implemented  in  Fortran  77  for  the  Silicon 
Graphics  (SGI)  workstation  architecture  and  given  the  code  name  ISOFIL. 
ISOFIL  was  developed  under  a  systems  integration  plan  based  oh  the  SGI 
Explorer  program  and  makes  use  of  extensions  to  Fortran  77,  including  the 
Fortran  POINTER  data  type  and  the  malloc  procedure  call.  ISOFIL  also 
includes  routines  from  two  public  domain  software  libraries,  BLAS  (level  1) 
[13]  and  SPARSPAK  [12].  All  performance  results  reported  in  this  section 
are  based  on  the  SGI  model  4D-35  workstation. 
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Numerical  validation  of  ISOFIL  so  far  has  included  evaluation  of  mass 
balance.  In  one  validation  exercise,  a  10  cm  radius,  center-gated,  disk  mold 
was  discretized  with  800  triangles  and  injected  at  the  constant  flow  rate 
of  one  cc/second  assuming  a  void  fraction  of  70%,  anisotropic  permeability 
(^11  =  1,A:22  =  .3),  and  two  fiber  orientations.  The  simulated  filling  time 
was  213.8  seconds  compared  to  the  expected  (analytical)  219.6  seconds,  or 
about  3%  relative  error.  The  same  model  was  evaluated  using  a  constant 
pressure  injection  of  one  kg/cm^.  The  mass  influx  was  approximated  by 
the  flux  across  the  control  volume  boundaries  surrounding  the  injector.  The 
resulting  approximation  was  214.1  cc  to  fill  the  mold,  versus  the  expected 
(analytical)  219.6  cc.  The  CPU  time  was  5  seconds,  including  input  and 
output. 

A  qualitative  evaluation  of  the  simulated  flow  fronts  indicates  that  flow 
front  shape  is  determined  by  element  shape.  This  result  can  also  be  inferred 
from  inspection  of  Figure  (2).  Highly  elongated  elements  lead  to  elongated 
control  volumes.  Elongated  volumes  may  fill  before  neighboring  volumes 
begin  to  fill,  in  a  locally  non-physical  manner.  This  effect  demonstrates  the 
need  for  a  relatively  uniform  mesh  with  a  good  aspect  ratio  for  the  triangles. 

The  performance  of  ISOFIL  was  compared  to  a  predecessor  code,  LIMS, 
from  the  University  of  Delaware  [7].  ISOFIL  uses  the  same  CVFE  approach 
as  LIMS  to  model  pressure  and  fluid  velocity.  There  are  minor  differences  in 
how  control  volume  flux  is  calculated.  The  primary  differences  are  that  LIMS 
assembles  and  factors  the  pressure  stiffness  matrix  and  the  flux  stiffness  ma¬ 
trix  at  each  time  step  of  the  CVFE  algorithm.  As  a  result,  the  performance 
differential  between  the  two  codes  increases  with  n,  the  number  of  nodes.  For 
a  2213-node,  4443-triangle  model  of  the  Ford  Aerostar  Crossmember  (see  be¬ 
low),  the  CPU  times  were  230  seconds  for  ISOFIL  and  22389  seconds  for 
LIMS  2.2,  a  speedup  factor  of  approximately  100. 

A  set  of  four  test  problems  are  presented  in  Figure  4,  including  a  plaque, 
disk,  auto  crossmember,  and  aircraft  keel  prototype.  The  plaque  and  disk 
models  are  simple  two-dimensional  geometries,  while  the  crossmember  and 
prototype  keel  box  are  thin-shell,  three-dimensional  structures.  These  mod¬ 
els  have  been  simulated  under  various  conditions,  including  different  choices 
of  injector/vent  locations,  material  type  (permeability  and  fiber  orientation). 
The  performance  results  are  presented  in  Table  1.  The  larger  test  problems 
require  several  minutes  of  CPU  time.  The  time  to  read  the  input  deck  (el¬ 
ement  mesh  and  material  properties)  is  included  in  the  results,  however  it 
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is  not  a  major  fraction  of  the  total  time.  The  error  in  mass  balance  ranges 
up  to  3%  for  the  test  problems  and  is  calculated  as  described  above  for  the 
constant  pressure  injection  case. 


Problem 

Nodes 

Elements 

Sec. 

Mass  Bal. 

Speedup 

disk 

442 

800 

5 

2.60% 

23 

plaque 

925 

1728 

22 

0.02% 

33 

xmbr 

2213 

4443 

230 

0.45% 

97 

proto 

2066 

4116 

188 

0.39% 

na 

Table  1:  ISOFIL  Performance  on  Four  Models.  Results  were  obtained 
on  an  SGI  4D-35,  CPU  time  includes  I/O.  Speedup  is  the  ratio  of  LIMS  2.2 
CPU  time  to  ISOFIL  CPU  time. 


4  Conclusions 

The  isothermal  RTM  filling  algorithm  has  a  computational  complexity  of 
(9(n^-®),  compared  to  more  general  RTM  algorithms  which  also  employ  the 
CVFE  formulation  but  require  O(n^)  computation.  This  relative  advantage 
results  in  a  100-fold  performance  improvement  over  a  similar  code  for  a  2213- 
node  model,  while  obtaining  the  same  solution  accuracy. 

The  speed  of  the  filling  simulation  is  critical  to  applying  simulation  results 
to  actual  mold  design.  A  fast  and  flexible  simulation  tool  allows  engineers  to 
include  modeling  in  the  design  process.  The  ISOFIL  code  is  currently  being 
used  for  interactive  filling  simulation  of  structural  aircraft  components  in 
connection  with  Army  procurement  projects.  The  complete  software  system 
permits  interactive  graphical  manipulation  of  mesh  and  material  properties, 
as  well  as  the  location  and  specification  of  injection  pressure/displacement 
time  profiles.  The  details  of  this  complete  system  will  be  published  in  a 
future  paper. 

Interactive  simulation  (and  real-time  control)  is  now  feasible  for  small 
(n  1,000)  problems  on  high  performance  workstations  and  supercomput¬ 
ers.  However  the  computational  requirement  of  operations  still  pro- 
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hibits  interactive  simulation  of  refined  three-dimensional  models  involving 
n  >  10®  nodes  (massively  parallel  supercomputers  are  capable  of  about  10^° 
floating  point  operations  per  second).  We  are  interested  in  simulations  which 
require  a  matter  of  seconds  or,  at  most  a  few  minutes,  of  real  time.  This  mo¬ 
tivates  the  need  for  further  investigation  of  filling  algorithms  which  depart 
from  the  conventional  CVFE  strategy. 
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Figure  4:  Four  Test  Problems  with  Flow  Front  Histories.  Clockwise 
from  upper  left:  disk,  crossmember,  keel  prototype,  plaque. 
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